You can download the materials for this notebook from this link.
First load the usual libraries
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.1.0; sf_use_s2() is TRUE
library(tmap)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)
library(tmaptools)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
Next, we’ll load a couple of datasets. These are from the week 7 lab materials.
abb <- st_read("abb-2193.gpkg")
## Reading layer `abb-2193' from data source
## `/home/osullid3/Documents/teaching/Geog315/slides/zoomed-in-maps/abb-2193.gpkg'
## using driver `GPKG'
## Simple feature collection with 2125 features and 16 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1742685 ymin: 5420357 xmax: 1785682 ymax: 5492218
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
welly <- st_read("wellington-base-data.gpkg")
## Reading layer `wellington-base-data' from data source
## `/home/osullid3/Documents/teaching/Geog315/slides/zoomed-in-maps/wellington-base-data.gpkg'
## using driver `GPKG'
## Simple feature collection with 122 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1735093 ymin: 5410973 xmax: 1774833 ymax: 5444626
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
These layers are in the same projection, so we don’t need to worry
about that, but as is usual, it is a good idea to transform all layers
to a matching projecting suitable for the project. Use
st_transform() to do this.
The easiest way to handle variations in scale and extent of layers is
a web map. For this, use tmap_mode("view")
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(welly) +
tm_polygons(col = "pop", alpha = 0.5) +
tm_shape(abb) +
tm_dots()
The different extents of these two datasets might be a problem. If we
are OK with filtering the point data to match the polygons data, then
just use st_filter():
welly_abb <- abb %>%
st_filter(welly)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(welly) +
tm_polygons(col = "pop", alpha = 0.5) +
tm_shape(welly_abb) +
tm_dots()
With a web map you can use the export options in RStudio to make snapshots of zoomed in portions of your maps for inclusion in a report.
Sometimes a web map is not a viable option. We have a number of static map options available. Before looking at them, I will switch back to static mapping:
tmap_mode("plot")
## tmap mode set to plotting
We can limit the extent of a static map with a bbox
option in the first layer. By default, tmap uses the extent
of the first layer to set the bounding box. So
tm_shape(welly) +
tm_polygons() +
tm_shape(abb) +
tm_dots()
makes a map that only extends as far as that first layer. If we want
it to extend to the full extent of the abb data then we can
do
tm_shape(welly, bbox = abb) +
tm_polygons() +
tm_shape(abb) +
tm_dots()
We’ll worry about the missing ‘base’ layer later in this document!
We can also use a bounding box to zoom in on an area. For example, if we make a bounding box that is from only the Lower Hutt City TA:
lower_hutt <- welly %>%
filter(ta_name == "Lower Hutt City")
tm_shape(welly, bbox = lower_hutt) +
tm_polygons() +
tm_shape(abb) +
tm_dots()
Again, there is missing basemap, but we’ll get to that in a bit.
If we want to make a specific bounding box, then we need to figure out the relevant coordinates. For this purpose, it is useful to make a temporary map with a grid.
tm_shape(welly) +
tm_polygons() +
tm_grid()
Now say we wanted a map showing only the central city. This has an
extent from 1,745,000 to 1,750,000 EW and 5,425,000 to 5,430,000 NS. We
can use this information to make a bounding box using the
bb function from tmaptools:
cbd <- bb(xlim = c(1745000, 1750000), ylim = c(5425000, 5430000))
and now we can do
tm_shape(welly, bbox = cbd) +
tm_polygons() +
tm_shape(abb) +
tm_dots()
There are other options in the bb function, that allow
you to specify it by using cx and cy (x and y
centre coordinates) and width and height. The
units are always the map projection units, which here are metres, so for
example another way to make the above map would be
cbd <- bb(cx = 1747500, cy = 5427500, width = 5000, height = 5000)
tm_shape(welly, bbox = cbd) +
tm_polygons() +
tm_shape(abb) +
tm_dots()
This approach is a bit easier to recalculate if you want to zoom in a
bit closer or zoom out a bit, by just changing the width and height of
the bounding box. There is also an ext option in the
bb function to increase or shrink the bounding box relative
to the requested width and height.
Often what we really want a web map for is the basemap.
We can get basemaps in static maps using the
rget_tiles() function in the package maptiles
and including it as a tm_rgb layer in the map. Here’s how
we use it. First we grab the map tiles:
library(maptiles) # you may have to install this
basemap <- get_tiles(abb)
Then we make a map
tm_shape(basemap) +
tm_rgb() +
tm_shape(abb) +
tm_dots()
If the basemap is a bit blurry, then you can zoom in a bit. You have
to be careful about this, as if you increase the zoom too much it will
take ages to download (and you won’t be able to read the detail anyway).
By experiment zoom levels in the range 9 to 11 are about right for the
Wellington region, and get_tiles is smart enough to get it
about right most of the time, without you specifying a zoom level. The
only other thing to watch out for is that maptiles delivers
whole web map tiles, so the cropping might be a bit strange. You can fix
this by setting the crop = TRUE parameter in the
get_tiles() function, or just using bbox in
tm_shape() as we have seen already:
# basemap <- get_tiles(abb, crop = TRUE) # this would crop the tiles to the data
tm_shape(basemap, bbox = abb) +
tm_rgb() +
tm_shape(abb) +
tm_dots()
If the colours are distracting, you can wash them out with the
saturation option in tm_rgb
tm_shape(basemap, bbox = abb) +
tm_rgb(saturation = 0.35) +
tm_shape(abb) +
tm_dots()
Other basemaps are available than. You can see the options in the
get_tiles() function help. For example, this is a nice one
(if a bit distracting):
basemap <- get_tiles(abb, provider = "Stamen.TerrainBackground")
tm_shape(basemap, bbox = abb) +
tm_rgb() +
tm_shape(abb) +
tm_dots()
Remember that if you are overlaying polygons on a basemap you will
probably want to set alpha (the opacity) to less than 1. So
here is a final example.
basemap <- get_tiles(welly)
tm_shape(basemap, bbox = welly) +
tm_rgb(saturation = 0.3) +
tm_shape(welly) +
tm_polygons(col = "pop", palette = "viridis", alpha = 0.5) +
tm_shape(abb) +
tm_dots() +
tm_layout(legend.outside = TRUE)
No really… here’s a final little extra. The cbd bounding
box might be something we want to use in a map like this. Unfortunately
a bb object is not able to be used directly to specify the
limits of a web map. We have to convert it to a simple features object
with a known projection. Here’s how we can do that.
cbd_sf <- cbd %>%
st_as_sfc() %>%
st_set_crs(2193) # this is the known projection in this case
basemap <- get_tiles(cbd_sf, zoom = 13) # I experimented to get best zoom
tm_shape(basemap, bbox = cbd) +
tm_rgb(saturation = 0.3) +
tm_shape(welly) +
tm_polygons(col = "pop", palette = "viridis", alpha = 0.5) +
tm_shape(abb) +
tm_dots() +
tm_layout(legend.outside = TRUE)